****************************************************************************
**		Authors: 	Or Tuttnauer & Liran Harsgor
**		Purpose: 	Produce Figure 4 in Kedar, Harsgor & Tuttnauer (JOP)  
**		input:		KHT_countrylevel.dta
*****************************************************************************

log using "figure 4.log", replace

set more off

use "KHT_countrylevel.dta", clear

* Generate x2,x3 with all values of modifier 
	set obs 902000
	gen magnitude_f = (int((_n-1)/200))/1100
	gen sr_f = (mod(_n-1,200))/36

** baseline model (Fig. 4a) **	
reg enps enpv lnmeddm vlmeddm smd mmm fused upper, robust

* Save coefficients and V-CV matrix *	
	matrix b=e(b)
	matrix V=e(V)
	scalar b1=b[1,1]  		/* The first coefficient (enpv) */
	scalar b3=b[1,3]		/* The first interaction's coefficient */
	scalar varb1=V[1,1]		/* The variance of b1 */
	scalar varb3=V[3,3]		/* The variance of b3 */
	scalar covb1b3=V[1,3]	/* The covariance of b1 and b3 */

* Marginal effect of enpv on enps
gen conbx_fict= b1+b3*magnitude_f
	
* Scatterplot *
recode conbx_fict (.00/.50=0) (.5001/.6=1) ///
		(.6001/.70=2)(.7001/.80=3)(.8001/.90=4)(.9001/1=5)(1.001/1.5=6), gen(conbx_cat)
		
	twoway  scatter sr_f magnitude_f if conbx_c==0 & sr_f>=1, mcolor(gs14) ///
			msize(vsmall) msymbol(D) || ///
			scatter sr_f magnitude_f if conbx_c==1 & sr_f>=1, mcolor(gs13) ///
			msize(vsmall) msymbol(D) || ///
			scatter	sr_f magnitude_f if conbx_c==2 & sr_f>=1, mcolor(gs12) ///
			msize(vsmall) msymbol(D) || ///
			scatter	sr_f magnitude_f if conbx_c==3 & sr_f>=1, mcolor(gs11) ///
			msize(vsmall) msymbol(D) || ///
			scatter	sr_f magnitude_f if conbx_c==4 & sr_f>=1, mcolor(gs9) ///
			msize(vsmall) msymbol(D) || ///
			scatter	sr_f magnitude_f if conbx_c==5 & sr_f>=1, mcolor(gs8) ///
			msize(vsmall) msymbol(D) || ///
			scatter	sr_f magnitude_f if conbx_c==6 & sr_f>=1, mcolor(gs6) ///
			msize(vsmall) msymbol(D) || ///
			scatter sr lnmeddm if dpr==1, msymbol(X) ///
			mcolor(black) mlabel(cyear) ///
			mlabcolor(black) mlabsize(vsmall) xtitle("Median magnitude (logged)") ///
			ytitle("Seat ratio") legend(off) ylabel(1 2 3 4 5)
	
graph save "effect_scatter_basesr.gph", replace
	
drop conbx_fict conbx_cat

** Our model (Fig. 4b) ** 
reg enps enpv lnmeddm sr vlmeddm vsr smd mmm fused upper, robust

* Save coefficients and V-CV matrix *	
	matrix b=e(b)
	matrix V=e(V)
	scalar b1=b[1,1]  		/* The first coefficient (enpv) */
	scalar b4=b[1,4]		/* The first interaction's coefficient */
	scalar b5=b[1,5]		/* The second interaction's coefficient */
	scalar varb1=V[1,1]		/* The variance of b1 */
	scalar varb4=V[4,4]		/* The variance of b4 */
	scalar varb5=V[5,5]		/* The variance of b5 */
	scalar covb1b4=V[1,4]	/* The covariance of b1 and b4 */
	scalar covb1b5=V[1,5]	/* The covariance of b1 and b5 */

* Marginal effect of enpv on enps
	gen conbx_fict= b1+b4*magnitude_f+b5*sr_f
	
* Scatterplot *
	recode conbx_fict (.00/.50=0) (.5001/.6=1) ///
		(.6001/.70=2)(.7001/.80=3)(.8001/.90=4)(.9001/1=5)(1.001/1.5=6), gen(conbx_cat)
		
	twoway  scatter sr_f magnitude_f if conbx_c==0 & sr_f>=1, mcolor(gs14) ///
			msize(vsmall) msymbol(D) || ///
			scatter sr_f magnitude_f if conbx_c==1 & sr_f>=1, mcolor(gs13) ///
			msize(vsmall) msymbol(D) || ///
			scatter	sr_f magnitude_f if conbx_c==2 & sr_f>=1, mcolor(gs12) ///
			msize(vsmall) msymbol(D) || ///
			scatter	sr_f magnitude_f if conbx_c==3 & sr_f>=1, mcolor(gs11) ///
			msize(vsmall) msymbol(D) || ///
			scatter	sr_f magnitude_f if conbx_c==4 & sr_f>=1, mcolor(gs9) ///
			msize(vsmall) msymbol(D) || ///
			scatter	sr_f magnitude_f if conbx_c==5 & sr_f>=1, mcolor(gs8) ///
			msize(vsmall) msymbol(D) || ///
			scatter	sr_f magnitude_f if conbx_c==6 & sr_f>=1, mcolor(gs6) ///
			msize(vsmall) msymbol(D) || ///
			scatter sr lnmeddm if dpr==1, msymbol(X) ///
			mcolor(black) mlabel(cyear) mlabcolor(black) mlabsize(vsmall) ///
			xtitle("Median magnitude (logged)") ylabel(1 2 3 4 5) legend(off)
	
graph save "effect_scatter_sr.gph", replace

graph combine effect_scatter_basesr.gph effect_scatter_sr.gph, commonscheme ///
ysize(10) xsize(19.5)scheme(s1mono) xcommon
graph save "figure 4.gph", replace

log close
